#!/usr/bin/env python
# -*- coding: iso-8859-2 -*-

import Image
import math
from ImageFilter import Kernel
from math import exp
from math import ceil
from math import sqrt
from __future__ import division

SMOOTH		= Kernel((3,3),[1,1,1, 1,1,1, 1,1,1])
SMOOTH_BIG	= Kernel((5,5),[1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1, 1,1,1,1,1])

SMOOTH_0		= Kernel((3,3),[1,1,1, 1,0,1, 1,1,1])
SMOOTH_BIG_0	= Kernel((5,5),[1,1,1,1,1, 1,1,1,1,1, 1,1,0,1,1, 1,1,1,1,1, 1,1,1,1,1])


LAPLAS_1_8	= Kernel((3,3),[-1,-1,-1,-1,8,-1, -1,-1,-1])
SOBEL_VER	= Kernel((3,3),[-1,0,1, -2,0,2, -1, 0, 1 ])
SOBEL_HOR	= Kernel((3,3),[-1,-2,-1, 0,0,0, 1, 2, 1 ])

BOOSTBLURFACTOR = 90.0
M_PI = 3.14159


def canny(image,sigma,tlow, thigh):
	imData = image.getdata()
	cols,rows = image.getdata()
	smoothedim = gaussian_smooth(imData, rows, cols, sigma)
	delta_x, delta_y = derrivativeXY(smoothedim, rows, cols)

	dir_radians = radianDirection(delta_x, delta_y, rows, cols, -1, -1) 
	magnitude = magnitudeXY(delta_x, delta_y, rows, col)


   non_max_supp(magnitude, delta_x, delta_y, rows, cols, nms);

   /****************************************************************************
   * Use hysteresis to mark the edge pixels.
   ****************************************************************************/
   if(VERBOSE) printf("Doing hysteresis thresholding.\n");
   if((*edge=(unsigned char *)calloc(rows*cols,sizeof(unsigned char))) ==NULL){
      fprintf(stderr, "Error allocating the edge image.\n");
      exit(1);
   }

   apply_hysteresis(magnitude, nms, rows, cols, tlow, thigh, *edge);

   /****************************************************************************
   * Free all of the memory that we allocated except for the edge image that
   * is still being used to store out result.
   ****************************************************************************/
   free(smoothedim);
   free(delta_x);
   free(delta_y);
   free(magnitude);
   free(nms);
}

def gaussianSmooth(image, rows,cols, sigma):
	kernel,windowsize = makeGaussianKernel(sigma)
	center = int(windowsize / 2)

   
   	#****************************************************************************
   	# Blur in the x - direction.
   	#***************************************************************************/
	for r in range(rows)::
		for c in range(cols):
         		dot = 0.0
         		sum = 0.0
         			for cc in range(-center), center+1):
            				if(((c+cc) >= 0) and ((c+cc) < cols)):
               					dot += (float)image[r*cols+(c+cc)] * kernel[center+cc]
               					sum += kernel[center+cc]
			tempim[r*cols+c] = dot/sum
      
   

	#****************************************************************************
	# Blur in the y - direction.
	#**************************************************************************/
	for c in range(cols):
		for r=0 in range(rows):
         		sum = 0.0;
         		dot = 0.0;
         		for(rr=(-center);rr<=center;rr++){
            			if(((r+rr) >= 0) && ((r+rr) < rows)){
               				dot += tempim[(r+rr)*cols+c] * kernel[center+rr];
               				sum += kernel[center+rr];
            
         
			smoothedim[r*cols+c] = int((dot*BOOSTBLURFACTOR/sum + 0.5))


def makeGaussianKernel(sigma):
	i = 0
	center = 0
	x = 0
	fx = 0
	sum = 0.0;

	windowsize = 1 + 2 * ceil(2.5 * sigma)
	center = int(windowsize) / 2)
	kernel = [None]*windowsize
	if DEBUG:
		print "windowsize: %s" % (windowsize, )


	for i in range (windowsize):
		x = i - center
		fx = pow(2.71828, -0.5*x*x/(sigma*sigma)) / (sigma * sqrt(6.2831853))
      		kernel[i] = fx
		sum += fx
  

	for i in range(windowsize):
		kernel[i] /= sum

	if DEBUG:
		print kernel
	
	return (kernel,windowsize)

def derrivativeXY(smoothedim, rows, cols):
	r = 0
	c = 0
	pos = 0

	for r=0 in range(rows):
		pos = r * cols
		delta_x[pos] = smoothedim[pos+1] - smoothedim[pos];
		pos = pos + 1
      		for c in range(1, cols-1):
         		(*delta_x)[pos] = smoothedim[pos+1] - smoothedim[pos-1]
			pos += 1
		delta_x[pos] = smoothedim[pos] - smoothedim[pos-1]

	for(c=0;c<cols;c++){
		pos = c;
		delta_y[pos] = smoothedim[pos+cols] - smoothedim[pos]
		pos += cols
		for r in range(1;rows-1):
         		delta_y[pos] = smoothedim[pos+cols] - smoothedim[pos-cols]
			pos += cols
		delta_y[pos] = smoothedim[pos] - smoothedim[pos-cols]
	return (delta_x, delta_y)

def radianDirection(delta_x,delta_y, int rows,int cols,int xdirtag, int ydirtag):
	r = 0
	c = 0
	pos = 0
	dirim = None
	dx = 0.0
	dy = 0.0
	
	dirim [None]*(rows*cols)
	dir_radians = dirim
	
	pos = 0
	for r in range(0,rows):
		for c in range(cols):
         		dx = (double)delta_x[pos];
         		dy = (double)delta_y[pos];

         		if (xdirtag == 1):
				dx = -dx
         		if (ydirtag == -1):
				dy = -dy
         		dirim[pos] = float(angle_radians(dx, dy))
			pos += 1
	

def angle_radians(x, y):
	xu = 0.0
	yu = 0.0
	ang = 0.0

	xu = abs(x)
	yu = abs(y)

	if((xu == 0) and (yu == 0)):
		return  0

	ang = atan(yu/xu);

   	if x >= 0:
      		if y >= 0:
			return ang
      		else:
			return 2*M_PI - ang
   	elif y >= 0:
		return M_PI - ang
      	else: 
		return M_PI + ang

def magnitudeXY(delta_x,delta_y, rows, cols):
	r = 0
	c = 0
	pos = 0
	sq1 = 0 
	sq2 = 0
	magnitude = [None]*(rows. cols)

	pos = 0
   	for r in range(rows):
		for c in range(cols):
			sq1 = int(delta_x[pos]) * int(delta_x[pos])
			sq2 = int(delta_y[pos]) * int(delta_y[pos])
			magnitude[pos] = int((0.5 + sqrt(sq1 + sq2))
	 		pos += 1
	return magnitude


non_max_supp(mag, gradx, grady, nrows, ncols) 
	rowcount = 0
	colcount = 0
	count = 0
	short *magrowptr,*magptr;
    	short *gxrowptr,*gxptr;
	short *gyrowptr,*gyptr,z1,z2;
	short m00,gx,gy;
	mag1 = 0.0 
	mag2 = 0.0
	xperp = 0.0
	yperp = 0.0
	unsigned char *resultrowptr, *resultptr;
    
	resultrowptr = result
	resultptr = result+ncols*(nrows-1)

   /****************************************************************************
   * Suppress non-maximum points.
   ****************************************************************************/
   for(rowcount=1,magrowptr=mag+ncols+1,gxrowptr=gradx+ncols+1,
      gyrowptr=grady+ncols+1,resultrowptr=result+ncols+1;
      rowcount<nrows-2; 
      rowcount++,magrowptr+=ncols,gyrowptr+=ncols,gxrowptr+=ncols,
      resultrowptr+=ncols){   
      for(colcount=1,magptr=magrowptr,gxptr=gxrowptr,gyptr=gyrowptr,
         resultptr=resultrowptr;colcount<ncols-2; 
         colcount++,magptr++,gxptr++,gyptr++,resultptr++){   
         m00 = *magptr;
         if(m00 == 0){
            *resultptr = (unsigned char) NOEDGE;
         }
         else{
            xperp = -(gx = *gxptr)/((float)m00);
            yperp = (gy = *gyptr)/((float)m00);
         }

         if(gx >= 0){
            if(gy >= 0){
                    if (gx >= gy)
                    {  
                        /* 111 */
                        /* Left point */
                        z1 = *(magptr - 1);
                        z2 = *(magptr - ncols - 1);

                        mag1 = (m00 - z1)*xperp + (z2 - z1)*yperp;
                        
                        /* Right point */
                        z1 = *(magptr + 1);
                        z2 = *(magptr + ncols + 1);

                        mag2 = (m00 - z1)*xperp + (z2 - z1)*yperp;
                    }
                    else
                    {    
                        /* 110 */
                        /* Left point */
                        z1 = *(magptr - ncols);
                        z2 = *(magptr - ncols - 1);

                        mag1 = (z1 - z2)*xperp + (z1 - m00)*yperp;

                        /* Right point */
                        z1 = *(magptr + ncols);
                        z2 = *(magptr + ncols + 1);

                        mag2 = (z1 - z2)*xperp + (z1 - m00)*yperp; 
                    }
                }
                else
                {
                    if (gx >= -gy)
                    {
                        /* 101 */
                        /* Left point */
                        z1 = *(magptr - 1);
                        z2 = *(magptr + ncols - 1);

                        mag1 = (m00 - z1)*xperp + (z1 - z2)*yperp;
            
                        /* Right point */
                        z1 = *(magptr + 1);
                        z2 = *(magptr - ncols + 1);

                        mag2 = (m00 - z1)*xperp + (z1 - z2)*yperp;
                    }
                    else
                    {    
                        /* 100 */
                        /* Left point */
                        z1 = *(magptr + ncols);
                        z2 = *(magptr + ncols - 1);

                        mag1 = (z1 - z2)*xperp + (m00 - z1)*yperp;

                        /* Right point */
                        z1 = *(magptr - ncols);
                        z2 = *(magptr - ncols + 1);

                        mag2 = (z1 - z2)*xperp  + (m00 - z1)*yperp; 
                    }
                }
            }
            else
            {
                if ((gy = *gyptr) >= 0)
                {
                    if (-gx >= gy)
                    {          
                        /* 011 */
                        /* Left point */
                        z1 = *(magptr + 1);
                        z2 = *(magptr - ncols + 1);

                        mag1 = (z1 - m00)*xperp + (z2 - z1)*yperp;

                        /* Right point */
                        z1 = *(magptr - 1);
                        z2 = *(magptr + ncols - 1);

                        mag2 = (z1 - m00)*xperp + (z2 - z1)*yperp;
                    }
                    else
                    {
                        /* 010 */
                        /* Left point */
                        z1 = *(magptr - ncols);
                        z2 = *(magptr - ncols + 1);

                        mag1 = (z2 - z1)*xperp + (z1 - m00)*yperp;

                        /* Right point */
                        z1 = *(magptr + ncols);
                        z2 = *(magptr + ncols - 1);

                        mag2 = (z2 - z1)*xperp + (z1 - m00)*yperp;
                    }
                }
                else
                {
                    if (-gx > -gy)
                    {
                        /* 001 */
                        /* Left point */
                        z1 = *(magptr + 1);
                        z2 = *(magptr + ncols + 1);

                        mag1 = (z1 - m00)*xperp + (z1 - z2)*yperp;

                        /* Right point */
                        z1 = *(magptr - 1);
                        z2 = *(magptr - ncols - 1);

                        mag2 = (z1 - m00)*xperp + (z1 - z2)*yperp;
                    }
                    else
                    {
                        /* 000 */
                        /* Left point */
                        z1 = *(magptr + ncols);
                        z2 = *(magptr + ncols + 1);

                        mag1 = (z2 - z1)*xperp + (m00 - z1)*yperp;

                        /* Right point */
                        z1 = *(magptr - ncols);
                        z2 = *(magptr - ncols - 1);

                        mag2 = (z2 - z1)*xperp + (m00 - z1)*yperp;
                    }
                }
            } 

            /* Now determine if the current point is a maximum point */

            if ((mag1 > 0.0) || (mag2 > 0.0))
            {
                *resultptr = (unsigned char) NOEDGE;
            }
            else
            {    
                if (mag2 == 0.0)
                    *resultptr = (unsigned char) NOEDGE;
                else
                    *resultptr = (unsigned char) POSSIBLE_EDGE;
            }
        } 
    }
}



if __name__ == '__main__' :

	import psyco
	psyco.full()
	canny_edges(1.5, 0.05, 1)
